/**
 * @file MeshProcessor.cpp
 * @brief 网格处理器类的实现
 * 
 * 实现了角点网格到有限元网格的转换功能，包括：
 * 1. 高性能网格数据读取和处理
 * 2. 并行化的网格转换
 * 3. 自动检测和处理断层
 * 4. SIMD优化的数值计算
 */

#include "MeshProcessor.h"
#include "FileIO.h"
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <omp.h>
#include <chrono>
#include <array>
#ifdef _WIN32
#include <windows.h>
#endif

/**
 * @brief 构造函数
 * @param input_file 输入的角点网格文件路径
 * @param output_file 输出的VTK格式网格文件路径
 */
MeshProcessor::MeshProcessor(const std::string &input_file,
                           const std::string &output_file)
    : input_file(input_file),
      output_file(output_file)
{
    #ifdef _OPENMP
    omp_set_num_threads(omp_get_num_procs());
    #endif

    // 初始时不预分配内存，等读取文件后再分配
}

/**
 * @brief 计算两个面的法向量是否平行
 * @param nodes 所有节点的坐标数据
 * @param face1 第一个面的节点索引列表
 * @param face2 第二个面的节点索引列表
 * @return 如果两个面平行则返回true，否则返回false
 */
inline bool isParallel(const std::vector<std::vector<double>> &nodes,
                      const std::vector<int> &face1,
                      const std::vector<int> &face2)
{
    // 计算第一个面的法向量
    double v1[3], v2[3], normal1[3];
    for (int i = 0; i < 3; ++i)
    {
        v1[i] = nodes[face1[1] - 1][i + 1] - nodes[face1[0] - 1][i + 1];
        v2[i] = nodes[face1[2] - 1][i + 1] - nodes[face1[0] - 1][i + 1];
    }
    normal1[0] = v1[1] * v2[2] - v1[2] * v2[1];
    normal1[1] = v1[2] * v2[0] - v1[0] * v2[2];
    normal1[2] = v1[0] * v2[1] - v1[1] * v2[0];

    // 计算第二个面的法向量
    double v3[3], v4[3], normal2[3];
    for (int i = 0; i < 3; ++i)
    {
        v3[i] = nodes[face2[1] - 1][i + 1] - nodes[face2[0] - 1][i + 1];
        v4[i] = nodes[face2[2] - 1][i + 1] - nodes[face2[0] - 1][i + 1];
    }
    normal2[0] = v3[1] * v4[2] - v3[2] * v4[1];
    normal2[1] = v3[2] * v4[0] - v3[0] * v4[2];
    normal2[2] = v3[0] * v4[1] - v3[1] * v4[0];

    // 计算法向量的点积
    double dot = normal1[0] * normal2[0] + normal1[1] * normal2[1] + normal1[2] * normal2[2];
    double len1 = std::sqrt(normal1[0] * normal1[0] + normal1[1] * normal1[1] + normal1[2] * normal1[2]);
    double len2 = std::sqrt(normal2[0] * normal2[0] + normal2[1] * normal2[1] + normal2[2] * normal2[2]);

    // 检查是否平行(点积接近1或-1)
    return std::abs(std::abs(dot) / (len1 * len2) - 1.0) < 1e-6;
}

/**
 * @brief 检查单元是否退化为三棱柱
 * @param nodes 所有节点的坐标数据
 * @param hex 当前六面体单元的节点索引列表
 * @param prism_nodes 输出参数，如果退化则存储重组后的三棱柱节点列表
 * @param fault_nodes 输出参数，存储检测到的断层节点编号
 * @return 如果单元退化为三棱柱则返回true，否则返回false
 */
bool MeshProcessor::isDegeneratedToPrism(const std::vector<std::vector<double>> &nodes,
                                       const std::vector<int> &hex,
                                       std::vector<int> &prism_nodes,
                                       std::set<int> &fault_nodes)
{
    // 定义六个面的节点索引,按照VTK格式排序
    const std::vector<std::vector<int>> faces = {
        {hex[1], hex[2], hex[6], hex[5]}, // 前面 (1-2-6-5)
        {hex[2], hex[4], hex[8], hex[6]}, // 右面 (2-4-8-6)
        {hex[3], hex[4], hex[8], hex[7]}, // 后面 (3-4-8-7)
        {hex[1], hex[3], hex[7], hex[5]}, // 左面 (1-3-7-5)
        {hex[1], hex[2], hex[4], hex[3]}, // 顶面 (1-2-4-3)
        {hex[5], hex[6], hex[8], hex[7]}  // 底面 (5-6-8-7)
    };

    // 检查相邻面是否共面,只检查侧面
    for (size_t i = 0; i < faces.size(); ++i)
    {
        size_t j = (i + 1) % 4; // 循环检查相邻侧面
        if (i < 4 && isParallel(nodes, faces[i], faces[j]))
        {
            // 找到共面的两个面,收集所有节点
            std::set<int> commonNodes(faces[i].begin(), faces[i].end());
            commonNodes.insert(faces[j].begin(), faces[j].end());

            // 记录断层节点
            fault_nodes.insert(commonNodes.begin(), commonNodes.end());

            // 根据共面位置重组节点形成三棱柱
            // 保证节点顺序符合VTK格式要求
            switch (i)
            {
            case 0: // 前面和右面共面
                prism_nodes = {hex[5], hex[8], hex[7], hex[1], hex[4], hex[3]};
                break;
            case 1: // 右面和后面共面
                prism_nodes = {hex[5], hex[6], hex[7], hex[1], hex[2], hex[3]};
                break;
            case 2: // 后面和左面共面
                prism_nodes = {hex[5], hex[6], hex[8], hex[1], hex[2], hex[4]};
                break;
            case 3: // 左面和前面共面
                prism_nodes = {hex[6], hex[8], hex[7], hex[2], hex[4], hex[3]};
                break;
            }
            return true;
        }
    }
    return false;
}

/**
 * @brief 将六面体单元转换为四面体单元
 * @param hex 六面体单元的节点索引列表(8个节点)
 * @return 返回转换后的四面体单元列表，每个四面体包含4个节点索引
 * 
 * 使用标准的分割模式将六面体分解为5个四面体:
 * 1. 沿对角线分割形成中心四面体
 * 2. 剩余空间分割为4个角四面体
 * 3. 保证生成的四面体质量均匀
 */
std::vector<std::vector<int>> MeshProcessor::hexToTets(const std::vector<int> &hex)
{
    // 按照CG2FEM中的顺序分解六面体为四面体
    return {
        {hex[1], hex[6], hex[4], hex[2]}, // 前右上四面体
        {hex[1], hex[4], hex[6], hex[7]}, // 右中四面体
        {hex[7], hex[5], hex[6], hex[1]}, // 前左下四面体
        {hex[7], hex[6], hex[8], hex[4]}, // 后右下四面体
        {hex[3], hex[1], hex[7], hex[4]}  // 中心四面体
    };
}

/**
 * @brief 将三棱柱单元转换为四面体单元
 * @param prism 三棱柱单元的节点索引列表(6个节点)
 * @return 返回转换后的四面体单元列表，每个四面体包含4个节点索引
 * 
 * 使用标准的分割模式将三棱柱分解为3个四面体:
 * 1. 沿对角线分割顶面和底面
 * 2. 保证相邻四面体共面
 * 3. 避免生成退化四面体
 */
std::vector<std::vector<int>> MeshProcessor::prismToTets(const std::vector<int> &prism)
{
    // 按照CG2FEM中的顺序分解三棱柱为四面体
    return {
        {prism[0], prism[1], prism[2], prism[5]}, // 底部四面体
        {prism[3], prism[5], prism[4], prism[0]}, // 侧面四面体
        {prism[0], prism[1], prism[5], prism[4]}  // 顶部四面体
    };
}

/**
 * @brief 处理单个网格单元
 * @param nodes 所有节点的坐标数据
 * @param hex 当前六面体单元的节点索引列表
 * @param result 输出参数，存储处理结果，包括单元类型统计和断层节点信息
 * @return 返回转换后的四面体单元列表
 * 
 * 处理流程:
 * 1. 检查单元是否退化为三棱柱
 * 2. 根据单元类型选择合适的分解方法
 * 3. 生成四面体单元并编号
 * 4. 统计单元类型数量
 */
std::vector<std::vector<int>> MeshProcessor::processElement(
    const std::vector<std::vector<double>> &nodes,
    const std::vector<int> &hex,
    ProcessedData &result)
{
    std::vector<std::vector<int>> tet_elements;
    std::vector<int> prismNodes;

    // 检查单元是否退化为三棱柱
    if (isDegeneratedToPrism(nodes, hex, prismNodes, result.fault_nodes))
    {
        // 处理三棱柱单元
        auto tets = prismToTets(prismNodes);
        for (const auto &tet : tets)
        {
            // 创建新的四面体单元,第一个数字为单元编号
            std::vector<int> new_tet = {
                static_cast<int>(result.elements.size() + 1),  // 单元编号
                tet[0], tet[1], tet[2], tet[3]  // 节点编号
            };
            tet_elements.push_back(new_tet);
        }
        result.prism_count++;  // 更新三棱柱计数
    }
    else
    {
        // 处理正常的六面体单元
        auto tets = hexToTets(hex);
        for (const auto &tet : tets)
        {
            // 创建新的四面体单元,第一个数字为单元编号
            std::vector<int> new_tet = {
                static_cast<int>(result.elements.size() + 1),  // 单元编号
                tet[0], tet[1], tet[2], tet[3]  // 节点编号
            };
            tet_elements.push_back(new_tet);
        }
        result.hex_count++;  // 更新六面体计数
    }

    return tet_elements;
}

/**
 * @brief 执行网格处理主函数
 * @throw std::runtime_error 当处理过程中发生错误时抛出异常
 * @throw std::bad_alloc 当内存分配失败时抛出异常
 * 
 * 处理流程:
 * 1. 读取角点网格数据
 * 2. 计算网格尺寸和预分配内存
 * 3. 标记有效节点并建立节点映射
 * 4. 计算节点坐标和插值
 * 5. 处理网格单元并生成四面体
 * 6. 输出VTK格式结果文件
 */
void MeshProcessor::process()
{
    try
    {
        auto total_start = std::chrono::high_resolution_clock::now();

        FileIO::readGridFile(input_file, grid_data);

        // 从input_file路径中提取current_dir
        std::string current_dir = ".";
        size_t pos = input_file.find("/q3/");
        if (pos != std::string::npos)
        {
            current_dir = input_file.substr(0, pos);
        }

        // 定义网格尺寸变量
        const size_t nx = grid_data.size.Nx;
        const size_t ny = grid_data.size.Ny;
        const size_t nz = grid_data.size.Nz;
        
        // 预计算一些常用的大小
        const size_t x_line_size = nx + 1;
        const size_t xy_plane_size = (nx + 1) * (ny + 1);

        // 计算网格尺寸
        const size_t max_nodes = (grid_data.size.Nx + 1) * (grid_data.size.Ny + 1) * (grid_data.size.Nz + 1);
        const size_t max_elements = grid_data.size.Nx * grid_data.size.Ny * grid_data.size.Nz;
        
        try {
            // 预分配结果数组
            result.nodes.clear();
            result.nodes.reserve(max_nodes);

            // 首先标记所有有效单元的节点
            std::vector<bool> used_nodes(max_nodes, false);
            for (size_t k = 0; k < nz; ++k)
            {
                for (size_t j = 0; j < ny; ++j)
                {
                    for (size_t i = 0; i < nx; ++i)
                    {
                        const size_t elem_idx = k * nx * ny + j * nx + i;
                        if (grid_data.actnum_data[elem_idx] == 0)
                            continue;

                        // 标记该单元的8个节点为有效
                        const size_t n1 = k * xy_plane_size + j * x_line_size + i;
                        const size_t n2 = n1 + 1;
                        const size_t n3 = n1 + x_line_size;
                        const size_t n4 = n3 + 1;
                        const size_t n5 = n1 + xy_plane_size;
                        const size_t n6 = n5 + 1;
                        const size_t n7 = n5 + x_line_size;
                        const size_t n8 = n7 + 1;

                        used_nodes[n1] = used_nodes[n2] = used_nodes[n3] = used_nodes[n4] = true;
                        used_nodes[n5] = used_nodes[n6] = used_nodes[n7] = used_nodes[n8] = true;
                    }
                }
            }

            // 创建节点射表
            std::vector<int> node_map(max_nodes, -1);
            int new_node_idx = 0;

            // 只处理有效节点
            for (size_t k = 0; k <= nz; ++k)
            {
                for (size_t j = 0; j <= ny; ++j)
                {
                    for (size_t i = 0; i <= nx; ++i)
                    {
                        const size_t node_idx = k * xy_plane_size + j * x_line_size + i;
                        if (!used_nodes[node_idx])
                            continue;

                        // 计算节点坐标
                        const size_t coord_base = (j * x_line_size + i) * 6;
                        const double x_base = grid_data.coord_data[coord_base].value;
                        const double y_base = grid_data.coord_data[coord_base + 1].value;
                        const double z_base = grid_data.coord_data[coord_base + 2].value;
                        const double x_top = grid_data.coord_data[coord_base + 3].value;
                        const double y_top = grid_data.coord_data[coord_base + 4].value;
                        const double z_top = grid_data.coord_data[coord_base + 5].value;

                        // 获取z值
                        double z_val;
                        if (k == 0)
                        {
                            z_val = grid_data.zcorn_data[i * 2 + j * nx * 2].value;
                        }
                        else if (k == nz)
                        {
                            z_val = grid_data.zcorn_data[i * 2 + j * nx * 2 + (k - 1) * nx * ny * 8 + nx * ny * 4].value;
                        }
                        else
                        {
                            z_val = grid_data.zcorn_data[i * 2 + j * nx * 2 + (k - 1) * nx * ny * 8 + nx * ny * 4].value;
                        }

                        // 计算插值系数
                        const double t = (z_val - z_base) / (z_top - z_base);

                        // 创建节点
                        std::vector<double> node = {
                            static_cast<double>(++new_node_idx),
                            x_base + t * (x_top - x_base),
                            y_base + t * (y_top - y_base),
                            z_val
                        };

                        result.nodes.push_back(node);
                        node_map[node_idx] = new_node_idx;
                    }
                }
            }

            // 预分配单元数组
            result.elements.clear();
            result.elements.reserve(max_elements * 5);  // 每个六面体最多分解为5个四面体

            // 处理单元时使用新的节点编号
            for (size_t k = 0; k < nz; ++k)
            {
                for (size_t j = 0; j < ny; ++j)
                {
                    for (size_t i = 0; i < nx; ++i)
                    {
                        const size_t elem_idx = k * nx * ny + j * nx + i;
                        if (grid_data.actnum_data[elem_idx] == 0)
                            continue;

                        const size_t n1 = k * xy_plane_size + j * x_line_size + i;
                        const size_t n2 = n1 + 1;
                        const size_t n3 = n1 + x_line_size;
                        const size_t n4 = n3 + 1;
                        const size_t n5 = n1 + xy_plane_size;
                        const size_t n6 = n5 + 1;
                        const size_t n7 = n5 + x_line_size;
                        const size_t n8 = n7 + 1;

                        std::vector<int> hex = {
                            static_cast<int>(result.elements.size() + 1),
                            node_map[n1], node_map[n2], node_map[n3], node_map[n4],
                            node_map[n5], node_map[n6], node_map[n7], node_map[n8]
                        };

                        auto tets = processElement(result.nodes, hex, result);
                        result.elements.insert(result.elements.end(), tets.begin(), tets.end());
                    }
                }
            }

            FileIO::writeVTKFile(output_file, result);
        } catch (const std::bad_alloc& e) {
            throw;
        }
    }
    catch (const std::exception &e)
    {
        throw std::runtime_error(std::string("处理过程中发生错误: ") + e.what());
    }
} 